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Stochastic modeling of reaction-diffusion kinetics has emerged as a powerful theoretical 
' ' tool in the study of biochemical reaction networks. Two frequently employed models are 

the particle-tracking Smoluchowski framework and the on-lattice Reaction-Diffusion Master 
Equation (RDME) framework. As the mesh size goes from coarse to fine, the RDME 
• initially becomes more accurate. However, recent developments have shown that it will 

become increasingly inaccurate compared to the Smoluchowski model as the lattice spacing 
^ becomes very fine. In this paper we give a new, general and simple argument for why the 

G RDME breaks down. Our analysis reveals a hard limit on the voxel size for which no local 

I I RDME can agree with the Smoluchowski model. 



1 Introduction 

A prevalent view in molecular systems biology is that the noise in cellular reaction networks, 
arising intrinsically from low copy numbers of macromolecules, can have a substantial impact on 
function [7, 10 . Two frequently used models for simulating stochastic reaction-diffusion systems 
are the reaction-diffusion master equation (RDME) [151 E] ^ind the Smoluchowski model [18] , 
which we will refer to as the mesoscopic and microscopic models, respectively. In the RDME the 
computational domain is divided into voxels. Diffusion is simulated as discrete jumps between 
adjacent vertices in the mesh. Reactions can occur locally when reactants are present in the same 
voxel. The RDME is attractive from a computational perspective; it is the logical extension of 
spatially homogenous simulations based on the Gillespie algorithm p!2|, and keeps track of the 
location of molecules only up to the resolution of the mesh, hence allowing for coarse-graining. 
Publicly available software packages for mesoscopic simulations based on the RDME include 
mesoRD [13j, SmartCeh pLj and URDME [5 . 

On a finer modeling level, the Smoluchowski model treats diffusion and reaction in continuous 
space, with molecules explicitly represented as spheres with a certain interaction radius. As such, 
it is an example of a model commonly referred to as particle-tracking. Software based on the 
Smoluchowski model are Smoldyn [2l[3l[4], MCell [16] and Green's function reaction dynamics 
(GFRD) [20l [21]. Smoldyn and MCell both discretize time whereas GFRD is continuous in 
time, making it more accurate at fine scales. 

A well known property of the mesoscopic model is that it converges to the classical reaction- 
diffusion partial differential equation in the macroscopic limit. For a system approaching the 
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microscopic regime, it is tempting to think of the RDME as a better and better approximation 
to the Smoluchowski model for finer and finer mesh resolutions. This picture is misleading, as 
in fact, it has been shown that as the size of the voxels in an infinite 3D domain decreases, 
all bimolecular reactions are eventually lost in the mesoscopic model [14 . Recent work using 
GFRD has demonstrated that fast, microscopic rebinding events can substantially affect the 
macroscopic properties of a biochemical signal cascade when reactions are highly diffusion limited 
[19 . To accurately simulate such systems requires a fine spatial resolution. On these scales, the 
conventional RDME may be too inaccurate to capture even the qualitative behavior predicted 
by the microscopic model [9 . 

The fact that the conventional mesoscopic model becomes inaccurate as we approach the 
microscopic level is not surprising, as we are moving out of the domain of validity for which it 
was derived. However it can pose a real practical problem, as it is hard to know a priori if a 
simulation with the RDME will yield useful or misleading results. This is especially true for 
biochemical models with multiscale properties, which are frequently encountered in molecular 
biology. Simply resorting to e.g. GFRD whenever in doubt is currently not feasible in general 
due to the high computational cost for systems with many particles. A natural approach to 
remedy this problem is to try to extend the domain of validity of the RDME as the mesh size 
tends to zero. Isaacson [14 suggests that one way of doing this would be to let the association 
rate constants depend explicitly on the meshsize. Recently, approaches to make such corrections 
to the RDME have been proposed [9l|8]. 

In this paper we show, by a simple argument, that below a certain critical size of the mesh 
it will be impossible to make the RDME consistent with the Smoluchowski model by local 
modifications to the bimolecular rates. Our result is valid for finite domains in 2D and 3D, and 
for mesh-dependent mesoscopic association rates. For a simple model problem we derive the 
correct mesoscopic association rates in 2D and 3D and compute the critical size of the mesh, 
and show that it can be considerably larger than the interaction radii of the molecules. In this 
light we review and discuss recent work in which the RDME has been modified in different ways 
in order to better agree with the microscale model for very small voxel sizes. 



2 Background 

2.1 The Reaction-Diffusion Master Equation 

In the mesoscopic model, the computational domain is divided into non-overlapping voxels. 
The state of the system is the discrete number of molecules of each biochemical species in 
each voxel in the mesh. Inside voxels, the molecules react according to a pre-specified set of 
rules. In a small time interval dt^ a bimolecular reaction A -\- B ^ for example, occurs with 
probability kaabdt/V^ where ka is the mesoscopic rate constant for the reaction, a and b are 
the copy numbers of A and B in that voxel, and V is the volume of the voxel. Diffusion is 
modeled as jumps between adjacent voxels. For a Cartesian mesh with mesh spacing /i, the rate 
for a diffusive jump is given by D/h^^ where D is the diffusion constant. The time evolution 
of the whole system is described as a continuous-time discrete-space Markov process and the 
probability density function (PDF) of the system evolves according to the RDME [151 HI]- 
Solving for the PDF directly is impossible in general due to the extremely large statespace, but 
realizations of the process can be generated efficiently using kinetic Monte Carlo methodology 

0. 
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2.2 The Smoluchowski model 



In the Smoluchowski model two molecules A and B are assumed to move by Brownian mo- 
tion with diffusion constants Da and Db, and react with a certain probability at a distance 
determined by the sum of their reaction radii, pA and ps- For two molecules A and 5, the 
probability of a bimolecular reaction is governed by the Smoluchowski equation. Given that 
the molecules start at positions xa^ and xb^ at time to, the equation for the PDF p of the 
new relative position r = xa — xb (in a spherical coordinate system r = (r,^, (/))), is given by 
dfT = DAt with initial condition p (r, to|ro, to) = ^ (r — ro) and boundary conditions 

lim p(r,t|ro,to) = 0, 

|r|^oo 



krP (r,t|ro,to) 



where p = pA ^ pB^ D = Da -\- Db and kr is the microscopic association rate. It can be 
shown that a weighted mean of the positions given by R = y^DsJ^AXA + \/ Da/ DbXb will be 
normally distributed [20 , and by sampling a new r and R we obtain the new positions of the 
molecules at some time t. 

An efficient method for simulating systems of molecules is the Green's function reaction 
dynamics (GFRD) [20l [21] method. In GFRD a system of molecules is decomposed into pairs 
of molecules by choosing a time step such that any given molecule is unlikely to react with more 
than one other molecule during that time step. This reduces the full problem to solving the 
Smoluchowski equation for pairs of molecules. 



3 Results 

3.1 Breakdown of the mesoscopic model 

Recent work has demonstrated that the RDME breaks down in the limit of infinitesimal voxels. 
Isaacson [14] shows that the probability for the occurrence of bimolecular reactions vanishes with 
decreasing voxel size for molecules on a lattice in an infinite 3D domain. The study is restricted 
to the case where the mesoscopic reaction rates are not dependent on the size of the voxels. Here 
we present a new and intuitive way to understand the degeneration of the mesoscopic model in 
finite domains in 2D and 3D, and with mesoscopic reaction rates that may or may not depend 
explicitly on the mesh. Our analysis will also provide additional insight into why and when this 
breakdown occurs. 

To see why the RDME model cannot work for very small voxels, it is illustrative to consider 
the simple process of bimolecular association 

A^B^C. (1) 
In 3D, the diffusion-controlled mesoscopic reaction rate ka is often taken to be 

This expression will hereafter be referred to as the conventional mesoscopic rate constant. It is 
valid for large enough voxels, and in 2D no analogous expression is well-defined. 

Consider the reaction ([T]) and, in a Cartesian coordinate system, let a molecule of A be 
stationary in some voxel Va, i.e. Da = 0, and let B diffuse with diffusion constant in a square 
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or cube with side length L and with periodic boundary conditions. Let Tmeso be the average time 
until the two molecules react in the discretized, mesoscopic model, let rj denote the average time 
for a diffusive jump and let k^eso be the mesoscopic reaction rate. Once the B molecule is in the 
voxel Va, the two molecules will react with probability /Cmeso/(^meso + '^j)^ or the B molecule 
diffuses away one voxel. Thus the two molecules will be on average M = {kmeso-\-'^j) /kmeso times 
in the same voxel before they react. Let td denote the average time until the two molecules are 
in the same voxel for the first time. The average time for an event (reaction or diffusion) given 
that the two molecules are in the same voxel is Tg = l/(/cmeso + 'Tj)^ and thus the average time 
for the B molecule to diffuse away one voxel and then back is Tg + r^, where is the average 
time until the B molecule reaches Va given that it starts one voxel away. This will be repeated 
on average M — 1 times until the molecules react. We can now write Tmeso as 

Tmeso = T/) + (M - 1) + T^) + Tg 

= rD^ ^meso(l + ^steps) =• + T^- (3) 

where N^t^p^ — '^J^h- following Theorem was proven in [17]. 

Theorem 1 Assume that the molecule B has a uniformly distributed random starting position 
xb on the lattice, xb not equal to xa, and that the molecules can move to nearest neighbors only 
(as in the RDME). Then as N ^ oo, the following holds: 

A^steps - ^-'N\og{N) + 0.19517V, 7V,\,p, ^ N (2D) 
iVsteps - 1.51647V, TVi^^p^ ^ TV, (3D). 

where N^teps ^-^ the average number of steps until xb = xa for the first time, N^t^p^ — '^J'^h 
the average number of steps until xb = xa given that A and B start one voxel apart and N is 
the number of voxels in the domain. 

Corollary 1 Let td be the time until A and B are in the same voxel for the first time. Then 

^3 

TD ^ ash^O (3D). 

where h is the voxel size. 

This follows immediately from Theorem [l] and the fact that tb = ^^^^ N^teps^J^ ^ where rj = 
2dD/h^ and d is the dimension. 

From Corollary [l] and Equation ([3| we conclude that for a sufficiently small voxel size in the 
discrete space model, we will have t meso ^ 'Td ^ micro (whcrc T micro the average time until 
the molecules react in the microscopic model), for any choice of the mesoscopic rate constant, 
since > for all /Cmeso > 0- Eventually, as ^ 0, no bimolecular reactions will occur since 
molecules can only react when they are in the same voxel and td ^ oo^ hence Tmeso ^ co. 

Note that the reason for this effect is not that the diffusion process is inaccurately described 
at these length scales. The problem lies in the assumptions that molecules can react only after 
having diffused to the same voxel. 



3.2 No local correction to the association rates can make the RDME 
consistent with the microscopic model. 

Due to the potential computational advantages of using the RDME over the microscopic model, 
it is natural to try to correct the RDME to agree better with the Smoluchowski model for fine 
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lattice spacings. The most obvious approach is to try to adjust the mesoscopic association rate 
constants in the RDME, making them dependent on the discretization. This would preserve the 
local nature of the reactions, and hence the low computational cost of the conventional RDME. 
One immediate consequence of the analysis in the previous section, however, is that below a 
certain mesh size h* no such local correction can make the mesoscopic mean association rate 
between two molecules agree with the microscopic model. In fact, for a given domain and model, 
this happens precisely when tjj > Tmicro- This is illustrated in Fig. [l]for the case of a square 
with side length L = 250 p and a cubic domain with side length L = lOOp, with kr ^ oo. 




Grid size h (multiples of p) Grid size h (multiples of p) 

(a) 2D (b) 3D 

Figure 1: The expected time until the molecules are in the same voxel for the first time, td, is 
computed with the RDME on a structured grid on a square (a) and a cube (b) with reflective 
boundary conditions. To the right of the vertical line we have rjj < Tmicro- In that region we 
could correct the mesoscopic reaction rate so that the expected time until the molecules react 
matches the expected time that we obtain by simulating the system on the microscale. To 
the left of the vertical line we have td > Tmicro , and it is impossible to change the mesoscopic 
reaction rate such that the result obtained with the RDME will match the result obtained from 
the Smoluchowski model. On the square we find that h* ^ 5.1p and on the cube h* ^ np. 
In our simulations we have used a uniform, structured Cartesian mesh and p = 2 ■ 10~^m, 
D = lO-i^m^s-i (3D) and D = lO-^W^-^ (2D). 

If this critical mesh size always occurred for < p, the problem would be less of a practical 
issue since the conventional RDME model makes little sense for mesh sizes smaller than the 
sum of the interaction radii of the molecules. However, in our example /i* occurs for /i* « irp 
(3D) and h* ^ b.lp (2D). As long as td < Tmicro^ that is > /i*, it is theoretically possible 
to modify the association rate, i.e. derive a discretization-dependent reaction propensity q{h) 
that will give the same mean association time as the microscopic model. Clearly, to be able to 
give the correct binding times for as large a regime as possible, the modified reaction constant 
should have the property q ^ oo for td ^ Tmicro- 

While our analysis does not preclude the possibility to better match the mean association 
time by increasing D and thus decreasing td^ this would make the effective diffusion too fast 
and thus introduce another source of error. 

4 Discussion 

Fig. [2] shows a schematic representation of the RDME's behavior as a function of the meshsize. 
For h < p^ the RDME makes little sense and we cannot expect the model to work in this regime. 
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Conventional RDME gives useful results 



No local correction 
possible 



Local correction possible, but model 
and geometry dependent 



Figure 2: Schematic representation of the RDME's behavior as a function of h. For h < ^ 
no local correction to the conventional mesoscopic reaction rates exists for the simple problem 
of diffusion to a target. 



For hmin < h < /imax the conventional mesoscopic rate constants will work well, but for h < hmin 
the RDME will become increasingly inaccurate. For h* < h < h^i^ it is possible to derive mesh 
and problem dependent reaction rates that make the RDME agree better with the microscopic 
model. The precise locations of the critical values /imin, ^max and /i* are model, geometry and 
discretization dependent. 

In Isaacson expands the analytical solutions of the RDME and the Smoluchowski equa- 
tion for the simple model problem A + 5 ^ in a series and computes the three leading terms 
in h. He shows that the two first terms will converge to the same value as h tends to zero, but 
that the difference between the third term will diverge. There is an h that minimizes the differ- 
ence between the first terms of the expansion. This does not strictly prove that, for sufficiently 
small /i, the error between the solutions increases as h decreases, and the result is valid in an 
unbounded domain, but it still illustrates that for some h < hmin the reaction rates will need to 
be modified to make the mesoscopic model accurate. As we have shown here, however, one will 
eventually reach /i* and the difference between the models will inevitably increase. 

Recently, two different corrections to ka have been proposed [El |9]. In the 3D case with 
a cubic domain and a uniform, Cartesian discretization, Erban and Chapman [8j consider the 

ka rh ki 

simple model problem A-\- B — ^ 5, — > A. They derive a mesh-dependent rate expression by 
matching the true steady-state distribution (which can be obtained analytically for this simple 
problem) to the distribution obtained using a meshsize h. They arrive at the following expression 
(in 3D) 

They also find a critical mesh size hcrit = Poo/{^aD) under which no further correction can 
be made to /c^, where j3oo ~ 0.25272 is a unitless constant valid for L ^ h. hcrit makes the 
denominator in Q zero and hence g ^ oo as ^ hcrit- Below that critical value of /i, the 
mesoscopic association reaction is defined to (7 = oo, i.e. association occurs as soon as the 
molecules are in the same voxel. Thus q satisfies the basic requirements of our analysis: the 
existence of a critical meshsize and the correct limiting behavior as the meshsize tends to that 
critical value. Substituting ka for the conventional expression and taking /c^ ^ 00 we obtain 
hcrit ^ ^P' Note that the propensity function Q was obtained without any comparison to the 
microscopic model. 

ka 

Fange et al. pursue a similar idea in [9 . They study the problem A -\- B ^ C, and 

kd 

derive mesoscopic reaction rates such that the equilibration time of the system matches the 
equilibration time in the Smoluchowski model. They carry out this analysis in 2D and 3D on a 
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disk and a sphere and obtain 



p{h) = kr/(l + a(l - I3)(l - 0.58/3)) (3D) 
p{h) = kr/{l + aln[l + 0.544(1 - (3)/ (3]} (2D) 



where /3 = p/{p -\- p -\- £ is the radius of a disk with area h'^ in 2D and a sphere with volume 
in 3D, a = kr/{4:7rpD) in 3D and a = kr/{27rD) in 2D. These expressions do not predict a 

critical mesh size, but has the property p{h) k^ as h ^ in both 2D and 3D. 

Based on our analysis we easily obtain another correction in both 2D and 3D. Equation Q 

suggest that we choose /cmeso to be 

^meso = (1 + ^steps)/ ('^micro — 'Td) 

in order to have Tmeso = Tmicro- For h sufficiently small we get 



From these expressions we obtain h* ~ y^exp(0.195l7r/2 + 3/4)p ^ 5.1p (2D) (where Tmicro 
has been approximated by the analytical expressions for a disk derived in [9]) and h* ^ 7rp (3D) 
(with Tmicro approximated by {ka/L^)~^) in good agreement with the simulations in Fig. 1. 

The corrections obtained by Erban and Chapman do not coincide with the corrections ob- 
tained by Fange et al., illustrating how the corrections are dependent on the ansatz used to 
derive them. On the other hand, our corrections given by (15) agree well with Erban and Chap- 
mans in 3D and predict a similar h* . As can be seen in Fig. 3j our corrections match the mean 
association time well in 2D (a) and all corrections give better results than the conventional ex- 
pression ka in 3D (b). Interestingly, Fange et al. find experimentally for their example that they 
cannot match the Smoluchowski model in [9l Fig. 3] perfectly using the conventional RDME 
with the local corrections p{h) even for h ~ 5p. Instead they modify the lattice model to allow 
for reactions between molecules in immediate neighboring voxels. In doing so they match the 
models all the way down to h = 2p. Our analysis explains why the local corrections alone were 
not sufficient, and their results demonstrate the possibility of better agreement with the micro- 
scopic model by a generalization of the conventional RDME to allow for neighbor-interactions. 

In this paper we have studied the behavior of the RDME for a simple but illustrative model 
problem: irreversible bimolecular association in the perfectly absorbing limit. For a finite but 
large kr (if the system is less diffusion-limited), /i* would shift to the left towards p in Fig. [l] 
However, for a reversible reaction microscopic reversibility must hold and if kr is finite then 
the mesoscopic dissociation rate will be larger than the microscopic dissociation rate for some 
h > . The /i* that we obtain here in the irreversible case will thus be a lower bound for 
the reversible case. In conclusion, the conventional RDME cannot be made consistent with the 
Smoluschowski model since there will always be a meshsize for which no local correction to the 
reaction rate can give the correct mean association time. Above /i* local corrections can extend 
the domain where the RDME works well. However, the corrections will inevitably be model and 
geometry dependent. 
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Grid size h (multiples of p) Grid size h (multiples of p) 

(a) 2D (b) 3D 

Figure 3: The expressions for the discretization-dependent mesocopic reaction rates from [SI [9] 
and those obtained here ah depend on the ansatz used to derive them. For irreversible association 
in 2D (a), our correction give a different association time than that of Fange et al. and in 3D 
(b) our correction agrees well with that of Erban and Chapman but gives a different result 
than that of Fange et. al. Despite this, all expressions produce more accurate results than the 
conventional expression for our model problem for h > . 
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